using DataFrames
using QuadGK
using Distributions
using JLD
using Optim
using ForwardDiff
using CSV
using LinearAlgebra
using SpecialFunctions
using Random
using DelimitedFiles
using SparseArrays
using Interpolations

clearconsole()

#-------------------------------------------------------------
# include Functions
#-------------------------------------------------------------
#cd("$(pwd())/Dropbox/Heterogeneity/Software/KS_Simulation/")
readDir = "$(pwd())/Functions/"
include(readDir *"logSpline_Procedures.jl");

#-------------------------------------------------------------
# Function for Gini coefficient
#-------------------------------------------------------------
function giniWiki(Xdata)
    Xsort  = sort(Xdata)
    Xfreq  = ones(size(Xsort))./length(Xsort)
    SX     = cumsum(Xsort.*Xfreq)
    GX     = SX[1]*Xfreq[1] + sum(Xfreq[2:end] .* (SX[2:end]+SX[1:end-1]))
    return 1 - GX/SX[end]
end

#-------------------------------------------------------------
# choose specification file
#-------------------------------------------------------------
nfVARSpec = "3"
specDir   = "$(pwd())/SpecFiles/"
include(specDir * "/fVARspec" * nfVARSpec * ".jl")

#-------------------------------------------------------------
# load data
#-------------------------------------------------------------
nDataSel = "3"
dataDir  = "$(pwd())/data/Para" * nDataSel *"/"

agg_data    = readdlm(dataDir * "aggvar_sim.csv", ',', Float64, '\n'; skipstart=0);
densdraws   = readdlm(dataDir * "a_cross_sim.csv", ',', Float64, '\n'; skipstart=0)';
employdraws = readdlm(dataDir * "empl_cross_sim.csv", ',', Float64, '\n'; skipstart=0)';
ss_K        = readdlm(dataDir * "ss_K.csv", ',', Float64, '\n'; skipstart=0);
K_exact     = agg_data[:,2] .+ ss_K;

# subsequently use the same knots regardless of sample size N,T
# knots_all = quantile(densdraws[employdraws.==1], quant_vec)

#-------------------------------------------------------------
# for each period compute quantiles
#-------------------------------------------------------------

vec_percs      = [0.1; 0.2; 0.5; 0.8; 0.9]
emp_percs_all  = zeros(T,1+length(vec_percs))
vec_cutoffs    = [1.0; 1.5; 2.0]
probMass_aboveX_all = zeros(T,1+length(vec_cutoffs))
gini_all       = zeros(T,2)
#timeidx = 1990+1/12 # 1990:M2

for tt = 1:T

    # time t data
    densdraws_t     = densdraws[1:N,tt]./K_exact[tt]
    employdraws_t   = employdraws[1:N,tt]
    selecteddraws_t = densdraws_t[employdraws_t.==1]

    emp_percs_all[tt,1]     = tt
    emp_percs_all[tt,2:end] = quantile(selecteddraws_t,vec_percs)

    probMass_aboveX_all[tt,1] = tt
    for ii = 1:length(vec_cutoffs)
        probMass_aboveX_all[tt,ii+1] = mean(selecteddraws_t .> vec_cutoffs[ii])
    end

    gini_all[tt,1] = tt
    gini_all[tt,2] = giniWiki(selecteddraws_t)

end

# check whether identical with the output from R (percentiles_data.csv)
CSV.write(dataDir * "percentiles_data_fVAR" * nfVARSpec * ".csv",DataFrame(emp_percs_all, :auto) )
CSV.write(dataDir * "probMass_aboveX_data_fVAR" * nfVARSpec * ".csv",DataFrame(probMass_aboveX_all, :auto) )
CSV.write(dataDir * "gini_coef_data_fVAR" * nfVARSpec * ".csv",DataFrame(gini_all, :auto) )
